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We describe the irrotational dust component of the universe in terms of a relativistic gradient 
expansion and transform the resulting synchronous metric to a Newtonian coordinate system. The 
two metrics are connected via a space-like displacement field and a time-like perturbation, pro- 
viding a relativistic generalization of the transformation from Lagrangian to Eulerian Coordinates. 
The relativistic part of the displacement field generates already at initial time a non-local density 
perturbation at second order. This is a purely relativistic effect since it originates from space-time 
mixing. We give two options, the passive and the active approach, on how to include the relativistic 
corrections for example in iV-body simulations. In the passive approach we treat the corrections 
as a non-Gaussian modification of the initial Gaussian field (primordial non-Gaussianity could be 
incorporated as well). The induced non-Gaussianity depends on scale and the redshift at which 
initial conditions are set, with /nl ~ few for small enough scales and redshifts. In the active ap- 
proach we show how to use the relativistic trajectory to obtain the initial displacement and velocity 
of particles for iV-body simulations without modifying the initial Gaussian field. 



I. INTRODUCTION 

Linearized cosmological perturbation theory (CPT) is 
a key technique for studying the nature of cosmological 
inhomogeneities [1, 2}. Its extension to second order has 
been applied to: (1) inflation, together with the subse- 
quent reheating [3-5], (2) to the coupled set of Einstein- 
Boltzmann equations of the primordial baryon-photon 
fluid [6-8], (3) to secondary effects after decoupling of 
the photons [9, 10], and (4) to late-time evolution of 
gravitational clustering [11, 12]. It is crucial to go to 
second order to understand non-linear aspects of the un- 
derlying physics, and to disentangle the various sources 
of non-Gaussianities. 

Here we report a somewhat different source of non- 
Gaussianity, relevant for the initial conditions of the 
CDM cosmic component, which cannot be directly em- 
bedded into one of the above groups (l)-(4). Instead, 
non-Gaussianity arises because of the non-linear coordi- 
nate transformation from a synchronous metric akin to 
Lagrangian coordinates to a Newtonian-like coordinate 
system, used for example in setting up JV-body simula- 
tions. In ref [13] it has been recently shown that such 
a coordinate transformation connects both metrics with 
a space-like displacement field and a time-like perturba- 
tion. The space-like displacement field is strongly related 
to the one in the Newtonian Lagrangian perturbation 
theory (LPT), and the time-like perturbation can be in- 
terpreted as the 3-velocity potential of the displacement 
field. This correspondence encourages us to interpret 
the synchronous coordinate system to be Lagrangian, 
attached to the CDM particles, whereas the observer's 
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"Eulerian" position is at rest in the Newtonian frame. In- 
deed, the metric perturbations in the Newtonian frame 
resolve to the Newtonian cosmological potential in the 
appropriate limit. 

Relativistic corrections in the displacement field should 
not influence the gravitational late-time evolution of the 
particle trajectories much. As we shall demonstrate how- 
ever, the relativistic coordinate transformation leads to 
non-Gaussian contributions in the density perturbations, 
which relate to the 3-velocity potential of the displace- 
ment field. These corrections matter already at the initial 
time and can be interpreted as a non-Gaussian modifi- 
cation of an otherwise initial Gaussian field. 

This paper is organized as follows. In section II we 
first present the gradient expansion metric in the syn- 
chronous frame for irrotational and pressureless cold dark 
matter (CDM) particles in a ACDM universe, where A 
denotes a cosmological constant. Then, we show how 
the residual freedom inherent in the gradient metric can 
be fixed, clarify its physical interpretation and separate 
the sources/origin of non-Gaussianities in the latter co- 
ordinate transformation, see section III. In section IV we 
interpret the relativistic corrections as the above men- 
tioned non-Gaussian modification of the initial Gaussian 
field ("the passive approach") and approximate them in 
terms of a bispectrum component; the particles are then 
displaced according to Newtonian LPT. In section V we 
describe the active approach where the initial Gaussian 
field is unaltered but the particles are displaced accord- 
ing to the relativistic trajectory. Then, in section VI we 
calculate the density contrast in the Newtonian coordi- 
nate system. Finally, we relate the density contrast in 
the Newtonian gauge to the one in iV-body simulations 
(section VII), and we conclude in section VIII. 
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II. THE GRADIENT EXPANSION METRIC 

We use the gradient expansion technique to solve the 
Einstein equations [14-19], although related relativistic 
approximation schemes lead to similar results [20-22]. 
We use a comoving/synchronous line element 



ds 2 = -At 2 + 7 W (t, q) Aq l Aq j , 



(1) 



where t is the proper time of the CDM particles and q 
are comoving/Lagrangian coordinates, constant for each 
pressureless and irrotational CDM fluid element. The 
gradient expansion approximates the spatial metric with 
a series containing an increasing number of spatial gradi- 
ents expressed in powers of the initial Riemann 3-tensor 
and its spatial derivatives. 1 Assuming standard inflation- 
ary initial conditions, the initial seed metric is 
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where <&(q) is the primordial Newtonian potential (in our 
case a Gaussian field), given at the initial time to- Using 
the formalism of ref [f 3] we then obtain up to four spatial 
gradients 2 
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where ", i" denotes a differentiation w.r.t. Lagrangian co- 
ordinate qi, summation over repeated indices is implied, 
and we have defined (A ^ 0): 
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with a(t) the scale factor. This result is valid for ACDM. 
To keep expressions as simple as possible we treat the 
Einstein-de Sitter (EdS) universe in the following, i.e., 

ft m = 1, 0\ = 0, and thus a(t) = {t/t a ) 2/3 . We shall 
generalize our findings to ACDM later. 

The precise limits in the above time-integrations mir- 
ror the chosen initial conditions which can be arbitrary. 
As we shall soon see the coefficients D and E give the 
time evolution of the displacement at linear and second- 
order respectively, and the velocity time coefficients are 
proportional to the time derivative of D and E. We thus 
need two constraints at any order, e.g., one for the initial 
displacement and the other for the initial velocity (al- 
ternatively one could replace one of them with the initial 
acceleration field) . To get the constraint for the displace- 
ment we require at initial time 
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This can be only achieved if D(to) = E(to) = 0, as can 
be easily verified through eq. (3). As we shall see the 
vanishing displacement field at initial time helps us to 
disentangle the sources of non-Gaussianity. 

The functions J and K, given by equations (5), con- 
strain the velocity coefficients proportional to D and E 
and depend on the actual physical situation (a dot de- 
notes a time derivative w.r.t. t). We give the general 
solution for D and E for generic initial conditions in ap- 
pendix A. Here we report a particular compelling class 
since their resulting expressions are closely related to the 
"slaved" initial conditions in refs [23-27] , and can thus be 
interpreted to be of the Zel'dovich type [28]. Since these 
slaved initial conditions require an initial velocity which 
is proportional to some acceleration, we think they model 
the nature of adiabatic fluctuations reasonably well. We 
require for this restricted class D(to) — E(to) = 0, 
I)(to) = 2t A and E(t ) = 2*g/21, which leads to 



D(t) = [a(t) - l]t 2 
3 



E(t) 



7 



a 2 (t) + a(t) 



(7) 



o ■ 



for eq. (3). The fastest growing modes in D and E do 
not change for any choice of initial conditions but the 
decaying modes do. Our conclusions do not change for 
any other realistic initial conditions and hence our find- 
ings do not depend on the specific choice of the decaying 
modes in (7); we make this choice to be concrete in what 
follows. Our final results are also valid for ACDM. 



III. NEWTONIAN COORDINATES AND THE 
DISPLACEMENT FIELD 



1 The Weyl tensor vanishes in three dimensions so the Riemann 
tensor is fully described in terms of the Ricci tensor. 

2 The same metric up to a transverse traceless tensor and without 
the decaying modes can also be found in equation (3.42) of [22]. 



To obtain a description in terms of the motion of par- 
ticles we transform the result of the gradient expansion 
from the comoving coordinates (*, q) to another coordi- 
nate system (t, x) . These two frames are connected by 
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the coordinate transformation 



according to (8), and 



Xi(t,q) = Qi +ipi(t,q) 
T(t,q)=t + C(t,q), 



(8) 



where tpi and C are supposed to be small perturbations. 
We shall soon identify xj) and £ as the displacement field 
and its velocity potential, respectively. We thus trans- 
form the metric 
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ds 2 = -dt 2 + ^(t, q) dq i dq j , 
to a Newtonian frame 



(9) 



ds 2 = -[l + 2A(T lX )]dT 2 

+ a 2 (r) [1 - 2B(t, x)} 5 ij dx i dx j , (10) 

where A, B are supposed to be small perturbations and 
7y is given in equation (3). Note that we have neglected 
the excitation of vector and tensor modes. This can 
be straightforwardly rectified if needed. The coordinate 
transformation requires: 
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We solve these equations for ip, C, A and B order by 
order. Formally, each small quantity is expanded in a 
series, i.e. A = eA^ + e 2 A^ + . . ., where e is supposed 
to be a small dimensionless parameter. The primordial 
potential $ is of order e. 

We report the results for the metric coefficients A and 
B in appendix B since they are not needed in the follow- 
ing. The coordinate transformation (8) — valid for arbi- 
trary initial conditions but restricted to an EdS universe, 
is up to second order 



Xi(t,q) 



2 



5Dd qi <£ 2 









01 







t2°2 : 



(14) 



where a dot denotes a time derivative w.r.t. t, the Hubble 
parameter for an EdS universe is H — 2/(3t), and we 
have defined 
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with 1/V 2 being the inverse Laplacian. 

Interestingly although not directly apparent, the ker- 
nel (18) has been also derived in ref [20] where the a of 
their eq. (18) is given by a = — 4V q 2 C 2 . In [20] a arises 
in the evolution equation of the second-order curvature 
perturbation in the Poisson gauge. 

The first line in eq. (14) agrees with Newtonian LPT 
[11, 25, 29, 30]; it contains the Zel'dovich approximation 
with its second-order improvement (2LPT). The rem- 
nant terms are relativistic corrections which should not 
bear much influence on the particle trajectories at late 
times. At initial time, however, they lead to an initial dis- 
placement and thus to a density perturbation. Also note 
that the time perturbations in the first line of eq. (15) cor- 
respond to the velocity perturbations in the Newtonian 
approximation. Indeed, v in the first line is the time co- 
efficient of the velocity potential at leading order in LPT, 
and the bracketed term in the same line in eq. (15) leads 
to oc g(t)V~ 2 G2 in an Eulerian coordinate system, with 
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which is the second-order kernel for the peculiar veloc- 
ity in Newtonian perturbation theory [11]. The remnant 
terms in (15) are again absent in Newtonian LPT. 

The precise choice of the decaying modes in D and E is 
of no great importance, reflecting only the chosen initial 
conditions [27]. However, it is very important to rec- 
ognize the disappearance of the Newtonian part of the 
displacement field in (14) for t — > i which occurs for 
any set of decaying modes. On the other hand, the last 
relativistic term in (14) is in general non-vanishing for 
t — > to. It generates the initial displacement (we addi- 
tionally take the divergence of the very equation) 



where we can easily read off the displacement field ip 
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For our simplified initial conditions, i.e., with the use of 
the growth functions (7) we have Vq — > t$. The time 
derivative of D is proportional to a velocity potential 
since D is the time coefficient of the displacement. Thus, 
the above expression vanishes only if the initial velocity 
between the synchronous and Newtonian frame vanishes; 
this initial data however would be unphysical since we 
already specified a vanishing displacement field. On the 
contrary, for a non-zero initial velocity the displacement 
kernel (18) is enhanced. 

Note again that the only restriction we have made to 
derive eq. (20) is to require jij(to,q) = fc„ initially. In 
a Newtonian treatment, i.e., a coordinate transformation 
which relates two Euclidean metrics, the requirement of 
7^ uclldean (io, q) — hj would imply an exact overlapping 
of the Eulerian and Lagrangian frames at initial time. 
Expression (20) would thus vanish. 

The above result is a purely relativistic effect. It results 
from the space-time mixing in the coordinate transfor- 
mation (11). Note that eq. (20) is also valid for ACDM 
with the corresponding D{t). 



IV. OPTION A: INITIAL NON— GAUSSIANITY 

The non-vanishing (relativistic) displacement field at 
initial time generates a density perturbation, which we 
shall use in the following to feed back to the Gaussian 
field: 
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with Ti = — (1 + 10ci/(3t^/ 3 )) and where c\ is a constant 
which depends on the chosen initial conditions (ci — > 
and so T\ — > — 1 for the above mentioned simplified initial 
conditions, see appendix A). The overall minus sign is the 
result of mass conservation. The last term in eq. (22) is 
the result of using the Newtonian time t and taking the 
limit t — ^ to for the relativistic trajectory (14) with the 
use of the time perturbation (15). 
In Fourier space eq. (22) becomes 
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where we have defined the (already symmetrised) non- 
local kernel 

FNL(fci,fc 2 ) =4v 2 Q [F^™(k 1 ,k 2 ) + F^(k ll k 2 )] , 

(24) 
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FIG. 1. Initial non-Gaussianity for the bispectrum (28) for 
different triangle configurations with fci23 = 0. We use the 
simplified initial conditions (with the respective growth func- 
tions (7)), and choose a redshift z — 50 such that cto — > 
5.5 Mpc/h. We plot the individual expressions in Fnl given in 
eq. (24), taking into account the non-Gaussian power induced 
through space-time mixing (STM) and through the time- 
perturbation (TP). The red (solid) curve corresponds to an 
equilateral triangle with k = ki = k2 = k^. The green (dashed) 
curve denotes a squeezed triangle with Afc = 0.012/i/Mpc and 
k = ki — &2 . The blue (dotted) curve denotes the contribu- 
tion through the time-perturbation and is identical for both 
squeezed and equilateral triangles. The yellow (shaded) re- 
gion specifies the total amount, STM+TP, of non-Gaussianity 
for squeezed configurations, and the black (double dotted) 
curve is the total amount of non-Gaussianity for equilateral 
triangles. 



where the kernels induced through space-time mixing 
(STM) and through the time perturbation (TP) are re- 
spectively 
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We use the short-hand notation k\ 2 = k\ + k 2 and k\ 2 = 
|fc 12 |. With — 2$ (see the following section) we have 



h = 5<P + <5, (2 



(27) 



and we recognize that the bispectrum induced through 
STM and TP is 

Bz{k 1 ,k2,k 3 ) = 2F NL {k 1 ,k2)P*(k 1 )P< s ,(k2) 

+ two permutations , (28) 
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at leading order, with the power spectrum 
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and similarily for the bispectrum. In the following we 
compare -Fnl from eq. (28) with the /^l parameter from 
the primordial non-Gaussianity of the local type. In fig. 1 
we plot the dimensionless quantity F NL . We have used 
the simplified set of initial conditions, with growth func- 
tions given in (7). We have specified the initial time to 
at a redshift z according to 
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for an EdS universe, where -fftoday = 100km s^/iMpc -1 , 
and h is the dimensionless Hubble parameter of order 
0.7. We have fixed z = 50 in the figure 3 but any other 
starting point is appropriate and just leads to a rescaling 
of the y-axis with the amplitude of -Fnl getting smaller 
with larger z. In this figure we show the results for equi- 
lateral and squeezed triangles (Afc = 0.012/i/Mpc). We 
have separated the non-Gaussian amplitude with respect 
to its origin, i.e., the amplitude from space-time mixing 
(STM) and from the time perturbation (TP). The yellow 
(shaded) region denotes the "sum" of both amplitudes for 
squeezed triangle configurations. Note that for squeezed 
triangles the amplitudes STM and TP are not partially 
summing/canceling itself as it is the case for equilateral 
triangles (which is denoted by the black [double dotted] 
curve). Instead, the kernel F^l (fci,fc 2 ) is vanishing for 
k\ = fc 2 and k\ + fc 2 + Afc = 0, and non-zero otherwise; 
this behaviour is reversed for the kernel F^ M . 

Equation (28) does only include the density perturba- 
tion induced through STM and TP (cf. eq. (22)). It is 
important to keep in mind that the initial acceleration 
between the frames (i.e., an initial acceleration of the 
particle) will induce a density perturbation as well — due 
to Einstein's equivalence principle. The appearance of 
an initial acceleration obviously depends on the chosen 
initial conditions. 

A precise scheme to implement this non-Gaussian 
component into an iV-body simulation can be found in 
ref [31, 32]. Then, the displacement and the velocity of 
the particles is just with respect to the Newtonian tra- 
jectory, i.e., it is given by Newtonian 2LPT [33]. 



The standard procedure to generate initial conditions 
in TV-body simulation is to use Newtonian 2LPT [33, 34]. 
In there, particles are usually displaced according to the 
fastest growing mode solutions of the second-order dis- 
placement field together with its velocity field. Whilst in- 
corporating relativistic corrections, however, one should 
include decaying modes as well, since they have roughly 
the same signature as the relativistic corrections. Addi- 
tionally and crucially, by neglecting the decaying modes 
one loses the control to adjust the initial velocities of the 
particles. As a consequence the magnitude of the initial 
velocity is rather accidental if decaying modes are ne- 
glected (since the initial velocity is then given by time 
differentiating of the fastest growing mode only). The 
resulting statistics of the matter fields are accidental as 
well, so artificial transiensts could flaw the (initial) grav- 
itational evolution. In appendix A we report the most 
general solution to adjust any initial data, depending on 
the actual physical situation. Here we describe how to 
set up generic initial data with the inclusion of (to be 
fixed) decaying modes together with the relativistic cor- 
rections. As before we require D(t ) = E{to) = 0. Note 
that this requirement does not harm the initial density 
configuration as long as the velocity and the acceleration 
field is generated according to the full displacement field. 

Suppose we have a simple rectangular grid in our N— 
body setting, thus without any spatial deformations. We 
adjust a universal clock which may be identified with an 
observers clock. We identify the universal clock to be 
r, given in eq. (15). Then, the particles position (14) is 
given by the initial grid position qn j k ^ = q plus the 
displacement at initial time (we suppress the i,j, k label 
for the specific grid points): 
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where c\ is a constant which depends on the chosen ini- 
tial conditions (see appendix A). This is the initial dis- 
placement of the particles from its original unperturbed 
grid positions. Note that at initial time the above can 
be either evaluated at the vicinity of the Eulerian or La- 
grangian spatial coordinate. This is however not true for 
the Eulerian and Lagrangian time. The velocity field at 
initial time is thus 



V. OPTION B: RELATIVISTIC TRAJECTORY 

Here we include the relativistic corrections in the so- 
called active approach. We call it active since not the 
Gaussian field is modified as in in the passive approach 
but the initial displacement of the particles. 



At z = 158 we have cto — > 1 Mpc/h. 
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where a prime denotes a differentiation with respect to 
Eulerian/Newtonian time t. 

Since we require a vanishing Newtonian displacement 
field at initial time we should specify the initial accelera- 
tion field for the sake to get a physical set of initial data. 
The obvious way to do so in an TV-body simulation is 
simply to require an initial acceleration component for 
each CDM particle. Instead of that one can take use of 
Einstein's equivalence principle and translate the accel- 
eration component into a density perturbation, i.e., allow 
different particle masses at different grid points. We re- 
port the required density perturbation in the following 
section. 

It is then straightforward to include the above into TV- 
body simulations. The general procedure how to obtain 
the displacement with its velocity on a grid can be found 
in appendix D2 in ref [33]. 



VI. THE DENSITY CONTRAST IN THE 
NEWTONIAN FRAME 

For the sake of completeness we also derive the density 
contrast in the new coordinate system. In the following 
it is convenient to restrict to the fastest growing modes 
only, since then we have a clean identification with the 
density contrast from the (Newtonian) literature. The 
spatial transformation of the fastest growing mode is 

X\(t, q)=q t + \a{t)tld^(t, q) + d^f (t, q) , (33) 



with the scalar 
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and the temporal transformation is 
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Note that the dependences in expressions (34) and (36) 
are either (t, q) or (r, x) at second order. The energy 
density p for an EdS universe written in the synchronous 
gauge is [13] 



p{t,q)=p{t) [l + 5(t,q)] 
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where p is the mean density. We thus transform the 
above according to (33) and (35). Then, we obtain the 



density contrast in the Newtonian coordinate system: 

S ^ X ^ = S^(t,x)+S^(t,x), (38) 
P 
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G2 is given in equation (19), and the Newtonian part of 
the densities at first and second-order are respectively 

/ \ 2 
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All dependences are w.r.t. (r, x), a vertical slash u de- 
notes a partial derivative w.r.t. Eulerian coordinate Xi, 
and we have neglected terms oc $ 2 which are not en- 
hanced by spatial gradients. It is also interesting to note 
the following relation: 
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(43) 



It links the second-order velocity perturbation to the 
second-order displacement field — through C 2 . Thus, a 
term oc at^Ci also arises in the density contrast (40). Its 
prefactor differs though in comparison with eq. (22) since 
we neglect decaying modes in this section but also since 
it is derived in a different gauge. 

Our result (38) seems to be in agreement with the 
second-order density contrast in the Poisson gauge in 
ref [20], i.e., we obtain the same spatial functions after 
a couple of manipulations. However, the prefactors seem 
to disagree. We leave the full analytic comparison of the 
second-order <5 in eq. (29) of ref [20] for a future project. 



VII. WHICH DENSITY IS MEASURED IN 
NEWTONIAN TV-BODY SIMULATIONS? 

In sections IV and V we have formulated two ways 
how to include the relativistic corrections. Both could 
be incorporated whilst setting up the initial conditions in 
TV-body simulations. In the last section we have calcu- 
lated the density contrast in the Newtonian gauge, so one 
might wonder which density will be measured in TV-body 
simulations. Since we assume that the onwarding grav- 
itational evolution in such simulations is still performed 
on a rectangular grid, the density contrast will be the 
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same as in the Newtonian approximation. If the simu- 
lation could entirely satisfy general relativity we would 
measure the density contrast in the Newtonian gauge. 
This is so since space-time is deformed in general rel- 
ativity and so is its 3-volume. Thus, to measure the 
(mass) density in the Newtonian gauge, eq. (38), we have 
to deform the volume Vg of the grid cells in the iV-body 
simulation according to 

Vg{t)= [ Jd 3 q, (44) 
Jg 

with the peculiar Jacobian [24, 35] 



J := sj detfffy]/ det[fe y ] , g iA := k ab x a A x b j , (45) 

where the second-order deformation tensor Xij is given 
by partial differentiation of eq. (14). The local density 
in the Newtonian gauge is then just the ratio of massive 
particles in the grid cell to the deformed volume Vg. 

VIII. CONCLUSIONS 

We obtained the relativistic Lagrangian displacement 
field together with its velocity potential from a gen- 
eral relativistic gradient expansion for an Einstein-de 
Sitter universe. By requiring that the metric reduces 
to a given initial condition in the synchronous comov- 
ing/Lagrangian frame we found a non-zero initial dis- 
placement field in the Newtonian/Eulerian frame which 
is absent in a purely Newtonian description. This initial 
displacement depends on the initial velocity potential of 
the particle, at second order in the primordial potential. 
We gave an example of initial conditions which are closely 
related to the so-called "slaved" ones and are thus of the 
Zel'dovich type. By including the decaying modes in the 
growth functions, we report the most general solution for 
the displacement and velocity in appendix A. 

The findings allow us to study the relevance of the ve- 
locity (of the fluid particle) between the synchronous and 
Newtonian coordinate system. We find that an initial ve- 
locity between the frames generates a non-local density 
perturbation. This is a purely relativistic effect which 
originates from two distinct sources, (1) space-time mix- 
ing in the coordinate transformation (11), and (2) from 
the fact that the particles' proper time is coordinate de- 
pendent. The net effect is not a result of the gravita- 
tional non-linear evolution, and it is even apparent if the 



initial acceleration is zero. (Non-zero accelerations gen- 
erate density perturbations even in the Newtonian treat- 
ment.) Equation (28) shows that the coordinate trans- 
formation induces a small velocity dependent amount of 
non-Gaussianity whose scale-dependence is depicted in 
fig. 1. Our result could be of importance in situations 
where the velocities of the traced objects are relatively 
high, since the /nl amplitude depends on the particle 
velocity potential squared. 

We considered the generation of initial conditions in 
./V-body simulations which we described in section IV: 
Instead of using the relativistic trajectory (14) one can 
treat the relativistic corrections in terms of an initial 
non-Gaussianity component, which acts as a correction 
to the (Lagrangian) potential and thus as a correction to 
the initial Gaussian field. We therefore adjust the cos- 
mological potential with respect to the initial non-local 
density perturbation but then let the particles evolve ac- 
cording to the Newtonian trajectory. This establishes a 
simple quasi-relativistic V-body simulation, and we call 
it the passive approach. In section V we described the 
active approach — the alternative way to incorporate the 
relativistic corrections in ./V-body simulations. In the ac- 
tive approach, the initial Gaussian field is unaltered but 
the particle trajectories include relativistic corrections. 
As an important note, we assume that the resulting ini- 
tial statistics in both approaches are equivalent, leave 
however the cumbersome proof for future work. If the 
resulting statistics agree, the same Lagrangian method 
could be used to generate non-Gaussian initial conditions 
for purely Newtonian ./V-body simulations. 

Finally, in section VII we showed how to relate the den- 
sity contrast from ./V-body simulations to the one in the 
Newtonian gauge. Essentially, space-time is deformed 
and so is the 3-volume of the grid cells. Thus, the den- 
sity in the Newtonian gauge is obtained by relating it to 
the appropriate physical volume. 
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Appendix A: General growth functions 

The general solutions for the displacement coeffi- 
cients (4) in an EdS universe are 



D(t) 
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where c\ - C4 are constants. To fix the constants one may 
choose specific initial data for D(to) and E(to) together 
with its first time derivatives. Again, this conforms to 
specify the initial displacement and the initial velocity. 

Generally it could be also appropriate to fix the initial 
peculiar gravitational field instead of the initial velocity, 
since the former sources a density perturbation and thus 
could flaw our main argument in this paper (i.e., that 
the initial non-Gaussianity is sourced only through ini- 
tial velocity perturbations). Indeed, defining the time 
coefficient of the initial peculiar gravitational field as 



ffpec(io) 



3i, 



D(t ) + D(t ), 



(A3) 



up to first order (similarity for the higher order coeffi- 
cients), we can construct displacements which are en- 
tirely without initial accelerations. We call this initial 
data "inertial", and the resulting displacement coeffi- 
cients are then 
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As can be proven by direct verification, the above modes 
do not change our conclusions: the coordinate transfor- 
mation yields to a density perturbation which is only 
sourced by a velocity perturbation. 



Appendix B: Perturbations in the Newtonian metric 

Here we summarise our findings for the metric 

ds 2 = -[l + 2A(T,x)]dT 2 

+ a 2 (r) [1 - 2B(t, as)] S tJ dx l dx 3 . (Bl) 

To avoid cumbersome expressions we restrict to an EdS 
universe and take only the fastest growing modes into 
account. Up to second order in $ the metric coefficients 
are 



A(t,x) 



(B2) 
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and 



with 



B(t,x) 



(B3) 



for this appendix. The term 0n can be obtained from a 
purely Newtonian treatment [17, 36], and it is just the 
cosmological potential up to second order. Indeed, mul- 
tiplying the second term on the RHS with V 2 /V 2 we 
find 



3 15 1 
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A vertical slash denotes a differentiation w.r.t. Eu- 
lerian coordinate Xi. Note that we have $ = &(to,x) 



^N = -<5>+la{T)t 2 -^F 2 (t ,x), (B7) 



with the well -known second-order kernel 
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F 2 (t ,x) = -$|«$| mm + <f>\l$\l m m + ^\lm$\lm- (B8) 



